Image enhancement

ABSTRACT

A method of deblurring images e.g. from an MRI scan, comprises calculating the distribution of the slope value of the intensities between the pixels in the image, and using a Levy distribution of power factor k of 1 (Cauchy distribution) up to but not including (2) to determine the correction term to be applied.

The present invention relates to a method for enhancing images, moreparticularly it relates to a method for the computer enhancement ofvisual images so that the original object can be more clearly seen andidentified.

Visual images which lack clarity are widely obtained for example fromCCTV cameras, from imaging of the body and brain e.g. mammography andother scans, from pictures obtained from telescopes e.g. in space orunder bad lighting conditions. A common problem with images is that thecontrast between two adjacent parts of the image is reduced or narrowedin the image compared with the original object and this means thatdetail is lost and the value of the image is reduced.

Another problem is the formation of false images such as “ringing”arising from an image enhancing process and/or pictures generated by animage processing method where artefacts such as noise or absent edgesoften deviate from the original.

A particularly important area is in the medical field where images ofthe human body or brain e.g. obtained by X-rays, ultra sound or byMagnetic Resonance Imaging (MRI) etc. have their value reduced by lackof detail.

In MRI of the brain there are several reasons why the images obtainedlack clarity and some of these e.g. movement of the patient, thelimiting resolution of the imaging or tomographic system and high noiselevels, which are encountered in functional MRI are likely to beinherent. It is important that the images can be enhanced and we havedevised a method for enhancing images.

In a known image processing method every image signal is replaced by amodified value which is based on the values of the image signals fromthe surrounding field of image elements. The signals from thesurrounding field are used to form a number of different combinationseach of which represents a different component of the image structurewithin the field.

Several methods are described in M. Bertero and P. Boccacci 1998Introduction To Inverse problems in Imaging (London: IOP).

These methods remove noise but introduce visible artefacts, particularlyin areas where the intensity values across the image are changingslowly, causing an artefact to appear as an edge.

We have now devised an improved method for enhancing images. This isbased on the fact that we have found that most images appear to obey aLevy distribution for their pixel slopes.

According to the invention there is provided a method for imageprocessing which comprises the acquisition or formation of an image,calculating the distribution of the slope value of the intensitiesbetween the pixels in the image, comparing this distribution with a Levydistribution of power factor k≧1, determining a correction term andapplying the correction term to the original image to obtain a processedimage.

Expressed as a Fourier transform, a zero centred, symmetric, Levydistribution contains two free parameters, viz scale and power (c,k) andis given by

${L_{c,k}(x)} = {\frac{1}{2\pi}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{- c}{\omega }^{k}}{\mathbb{e}}^{{\mathbb{i}\omega}x}{\mathbb{d}\omega}}}}$

For the Levy distributions used in the present invention the powerfactor (k) can range from 1 (Cauchy Distribution) up t⁻2 (Gaussiandistribution). In different images different power factors in the Levydistribution can be exploited, the power factor (k) chosen beingselected on a priori information elucidated by prior analysis and/orexperience to give the best fit.

The method can be repeated iteratively using the image result from aformer iteration as the starting point image for a subsequent iteration.A preferred realisation involves holding the non-fixed parameter c fromthe first iteration to completion. Other iterations allow the iterativeupdating of c.

The image can be formed from any conventional method in which an imageis obtained and converted to electronic signals, normally digitised, andthen it can be projected onto a screen e.g. a VDU or to a printer orother device for inspection. The method of the present invention enablesthe image projected onto the screen to be enhanced so as to present animage of greater clarity than otherwise be achieved.

The invention works for many types of image and is particularlyapplicable to images obtained from MRI.

A Cauchy distribution gives good results even for images with Levyindex >1, this is because significant enhancement may still be achievedeven when the distribution is not tailored to the image statistics,whereupon further enhancements are achieved. The calculation of theslopes and obtaining the distribution of the slopes can be carried outby a computer using known techniques.

When an image is blurred during capture by an imaging system of unknowncharacteristics, this invention, in conjunction with a blinddeconvolution method may both resolve the quality of the image andcharacterise the imaging properties of the image capture system. Theimage can be deconvolved by using known methods e.g. as described in“Blind Deconvolution by Means of the Richardson-Lucy Algorithm” D AFish, A M Brinicombe, E R Pike and J G Walker Journal of The OpticalSociety of America A, 12, 58-65 1995.

The invention is illustrated in detail to reconstruct an unknown objectthat has been blurred by a space invariant aberration. Using the Fouriermethod to model the problemg(i,j)=f(i,j)

k(i,j)+η(i,j)

Where, g(i,j) is the measured data, k(i,j) is the point spread function,ƒ(i,j) is the unknown object and η(i,j) is a noise process.

The slopes of the image can be approximated in the x and y directions bya finite difference scheme. Using the forward differenceδ_(x)(x,y)=ƒ(x+1,y)−ƒ(x,y)δ_(y)(x,y)=ƒ(x+1,y)−ƒ(x,y)

Setting the total slope, s as the magnitude of the combined δ_(x) andδ_(y) slopes gives

$\begin{matrix}{{s\left( {x,y} \right)} = \sqrt{\delta_{x}^{2} + \delta_{y}^{2}}} \\{= \sqrt{\left\lbrack {{f\left( {{x + 1},y} \right)} - {f\left( {x,y} \right)}} \right\rbrack^{2} + \left\lbrack {{f\left( {x,{y + 1}} \right)} - {f\left( {x,y} \right)}} \right\rbrack^{2}}}\end{matrix}$

For example using a two dimensional Levy distribution with k=1, known asthe symmetric Cauchy distribution in R², as the distribution to describethe probability density function of the (δ_(x),δ_(y)) slopes thecalculation can be progressed analytically. For other values of knumerical methods must be used. Thus for k=1

$\begin{matrix}{{P\left( {\delta_{x},\delta_{y}} \right)} = {\frac{1}{2\;\pi}\frac{c}{\left( {c^{2} + \delta_{x}^{2} + \delta_{y}^{2}} \right)^{3/2}}}} \\{= {\frac{1}{2\;\pi}\frac{c}{\left( {c^{2} + s^{2}} \right)^{3/2}}}}\end{matrix}$

Assuming that all the increments are uncorrelated the likelihood of theimage is a product of the probability density of all pixel values.Taking the logarithm of this value (the log likelihood).

$\begin{matrix}{{\log(L)} = {\log\left\{ {\prod\limits_{x,y}^{\;}\;{P\left( {{\delta_{x}\left( {x,y} \right)},{\delta_{y}\left( {x,y} \right)}} \right\}}} \right.}} \\{= {\log\left\{ {\prod\limits_{x,y}^{\;}{\frac{1}{2\;\pi}\frac{c}{\left( {c^{2} + {\delta_{x}^{2}\left( {x,y} \right)} + {\delta_{y}^{2}\left( {x,y} \right)}} \right)^{3/2}}}} \right\}}} \\{= {\sum\limits_{x,y}{\log\left\{ {\frac{1}{2\;\pi}\frac{c}{\left( {c^{2} + {s^{2}\left( {x,y} \right)}} \right)^{3/2}}} \right\}}}} \\{= {\sum\limits_{x,y}{{- \log}\left\{ {\frac{2\;\pi}{c}\left( {c^{2} + {s^{2}\left( {x,y} \right)}} \right)^{3/2}} \right\}}}}\end{matrix}$

It is desired to maximise Q whereQ=log(L)−λ×MSEand MSE is the mean square error given by

${MSE} = {\sum\limits_{i,j}\left\{ {{g\left( {i,j} \right)} - {{f\left( {i,j} \right)} \otimes {k\left( {i,j} \right)}}} \right\}^{2}}$

From this the value of λ can be determined which will maximise Q as

$\left. \Rightarrow{d\left( {i,j} \right)} \right. = {3\left\{ {\frac{{\delta_{x}\left( {i,j} \right)} + {\delta_{y}\left( {i,j} \right)}}{c^{2} + {s^{2}\left( {i,j} \right)}} - \frac{\delta_{x}\left( {{i - 1},j} \right)}{c^{2} + {s^{2}\left( {{i - 1},j} \right)}} - \frac{\delta_{y}\left( {i,{j - 1}} \right)}{c^{2} + {s^{2}\left( {i,{j - 1}} \right)}}} \right\}}$andD(u,v)+λ{G(u,v)−F(u,v)K(u,v)} K ⁻(u,v)=0

After solving for F(u,v) this becomes

${F\left( {u,\upsilon} \right)} = {\frac{{G\left( {u,\upsilon} \right)}{K^{*}\left( {u,\upsilon} \right)}}{{{K\left( {u,\upsilon} \right)}}^{2}} + {\frac{1}{\lambda}\frac{D\left( {u,\upsilon} \right)}{{{K\left( {u,\upsilon} \right)}}^{2}}}}$

Taking a solution

$\begin{matrix}{{B\left( {u,\upsilon} \right)} = \frac{D\left( {u,\upsilon} \right)}{{{K\left( {u,\upsilon} \right)}}^{2}}} \\{{f_{new}\left( {i,j} \right)} = {{f_{old}\left( {i,j} \right)} + {\frac{1}{\lambda}{b\left( {i,j} \right)}}}}\end{matrix}$the total slopes will be updated by

${s_{new}\left( {i,j} \right)} = {{s_{old}\left( {i,j} \right)} + {\frac{1}{\lambda}{s_{b}\left( {i,j} \right)}}}$where b(i,j) is a correction term so that λ is given by

$\lambda = {- \frac{\sum\limits_{i,j}\frac{s_{b}^{2}\left( {i,j} \right)}{c^{2} + {s_{new}^{2}\left( {i,j} \right)}}}{\sum\limits_{i,j}\frac{{s_{old}\left( {i,j} \right)}{s_{b}\left( {i,j} \right)}}{c^{2} + {s_{new}^{2}\left( {i,j} \right)}}}}$

A similar analysis can be performed using a centred difference formulato approximate the slope values which leads to a slightly differentexpression for d(i,j)

The process is preferably reiterated a number of times with the processrepeated on the previously corrected image to obtain furtherenhancement, preferably there is alternating use of forward differenceand centred difference in each iteration.

In one preferred embodiment, an image is made by some means and a Wienersolution is calculated. The total slope histogram for the reconstructionis calculated and fitted by a symmnetric Cauchy in R². From this acorrection term c is obtained.

Using the x,y and s slopes and c, the correction term written d(i,j) iscalculated. This term is transformed into Fourier space and used toobtain B(u,v). {circle around (2)} is then calculated and, afterdividing B(u,v), is added to the old reconstruction, the newreconstruction is then put back into real space and the processreiterated.

The number of iterations used to obtain the best results can bedetermined from a priori assumptions about the appearance of theoriginal object.

A similar analysis may be performed using a centred-difference formulato approximate the slope values.

Thus it is possible to enhance the original image.

The reconstruction formed does not have to have a slope which closelymatches a Levy distribution and it has been found that the method of theinvention works on synthetic shapes.

The method of the invention was demonstrated by taking an image of anoriginal object, blurring it by using a blurring function andreconstructing it using a conventional method (Wiener reconstruction)and by using the method of the invention, alternating forward andcentred difference, k=1; the results are illustrated in the accompanyingdrawings in which:—

FIG. 1 a shows the original object,

FIG. 1 b shows the blurring function

FIG. 1 c shows the recorded blurred and noisy image

FIG. 1 d shows a conventional reconstruction

FIG. 2 a shows a reconstruction using the method of the invention afterone iteration

FIG. 2 b after three iterations

FIG. 2 c after six iterations and

FIG. 2 d after 20 iterations.

As can be seen the method of the invention gives a superiorreconstruction of the object.

The invention claimed is:
 1. A method for image processing which comprises forming an image, calculating the distribution of the slope value of the intensities between the pixels in the image, comparing this distribution with a Levy distribution of power factor k≧1, determining a correction term and applying the correction term to the image to obtain a processed image.
 2. A method as claimed in claim 1 in which the power factor is from 1 up to, but not including,
 2. 3. A method as claimed in claim 1 in which the power factor is
 1. 4. A method as claimed in claim 1 in which the method is repeated iteratively using the processed image from a prior iteration as the starting image for a subsequent iteration.
 5. A method as claimed in claim 4 in which a non-fixed parameter c defining the distribution for the first iteration is used in subsequent iterations.
 6. A method for image processing which comprises forming an image, subjecting the image to blind deconvolution to obtain a blurring function, calculating the distribution of the slope value of the intensities between the pixels in the image, comparing this distribution with a Levy distribution of power factor k≧1, determining a correction term and applying the correction term to the image to obtain a processed image.
 7. A method as claimed in claim 1 in which the image is obtained is an image of part or all of the human body.
 8. A method as claimed in claim 7 in which the image is obtained by X-rays, ultra sound or by magnetic resonance imaging.
 9. A method as claimed in claim 6 in which the power factor is from 1 up to, but not including,
 2. 10. A method as claimed in claim 6 in which the power factor is
 1. 11. A method as claimed in claim 6 in which the method is repeated iteratively using the processed image from a prior iteration as the starting image for a subsequent iteration.
 12. A method as claimed in claim 11 in which a non-fixed parameter c defining the distribution for the first iteration is used in subsequent iterations. 